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Abstract 

In the e-regime of QCD the main features of the spectrum of the low-lying eigenvalues of 
the (euclidean) Dirac operator are expected to be described by a certain universality class 
of random matrix models. In particular, the latter predict the joint statistical distribution 
of the individual eigenvalues in any topological sector of the theory. We compare some of 
these predictions with high-precision numerical data obtained from lattice QCD for a range 
of lattice spacings and volumes. While no complete matching is observed, the results agree 
with theoretical expectations at volumes larger than about 5fm*. 



1. Introduction 

The proposition that the low-lying eigenvalues of the Dirac operator in the so-called 
e-regime of QCD [1,2] are distributed in the same way as the eigenvalues of a large 
random matrix was put forward a number of years ago [3-5] and has since then 
been worked out in great detail (see refs. [6,7] for a review and further references). 
Perhaps the most important qualitative prediction of random matrix theory is that, 
at vanishing quark masses, the eigenvalues scale proportionally to where S 

denotes the (bare) quark condensate and V the space-time volume. In particular, 
the spectrum near the origin rapidly becomes very dense when the volume increases. 

* On leave from CNRS, Centre de Physique Theorique, F-13288 Marseille, France 
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So far the correctness of the proposition has not been estabhshed from first prin- 
ciples. Apart from symmetry and universahty arguments, it is supported by chiral 
perturbation theory, where the Lcutwyler-Smilga sum rules [2] can be shown to be 
reproduced, at large volumes, by random matrix theory. To date there is, however, 
no similar theoretical check on the distributions of the individual eigenvalues [8,9]. 

Random matrix theory is well defined for any number Ngea > of sea quarks. It 
is thus possible to extend the proposition to quenched QCD, even though the latter 
tends to be singular in the chiral limit. Another instance where random matrix 
theory with iVgea = may be expected to apply is full QCD with quark masses 
m such that mHV S> 1. The quark determinant is safely bounded from below in 
this case, and the distributions of the low-lying eigenvalues of the massless Dirac 
operator should consequently be as in quenched QCD (up to a scale transformation). 

The present paper is part of an ongoing project whose final goal is to extract 
physical parameters, such as the pion decay constant and the electroweak effective 
couplings, from numerical simulations of lattice QCD in the e-regime [10,11]. Since 
chiral symmetry plays a central role in this context, we use a lattice Dirac operator 
that satisfies the Ginsparg-Wilson relation and thus preserves the symmetry [12-20]. 
Unfortunately numerical simulations of full QCD then become even more demanding 
than they normally are, and in this paper we shall, therefore, only consider the case 
of quenched QCD. 

Comparisons of random matrix theory with simulation data obtained from lattice 

QCD with exact chiral symmetry have previously been published by a number of 
collaborations [21-24]. Here we extend these studies to significantly larger lattices 
and collect data at different lattice spacings so as to be able to check for lattice 
effects. We wish to add that simulations in the e-regime of QCD become increasingly 
difficult at large volumes, because the low- lying eigenvalues of the Dirac operator 
are extremely small and closely spaced. The use of efficient techniques such as those 
described in ref. [10] is thus essential. 



2. Random matrix model 

To set up notations, we now briefly recall the definition of the matrix model that is 
expected to describe the spectrum of the Dirac operator in QCD. There is actually 
a whole universality class of such models, and we will only describe the simplest 
representative, the so-called gaussian chiral unitary model. 
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Let us consider N x N matrices of the form 



where is a complex rectangular random matrix. We use the symbol D to indicate 
that the matrix represents the massless Dirac operator in the matrix model. It is 
then natural to think of N as the space-time volume (times some proportionality 
constant), while the block structure of D is interpreted as a chiral decomposition. 
Moreover, since any matrix of this form has \N^ — N- \ chiral zero modes, the index 

iy = N+-N_ (2.2) 

may be identified with the topological charge in QCD. 

The observables O in the matrix model are arbitrary functions of matrix elements 
of W, and we are interested in the large- A'^ limit of the expectation values 

(O), = ^ / ^W] O det(D + m)^- g-|jvtr{wtw} ^2.3) 
J 

at fixed /X = mN. Clearly m represents the quark mass in the matrix model (in the 
case where all quarks have the same mass), but this parameter is actually irrelevant 
in the present paper, since we shall compare the matrix model with quenched QCD 
and thus set A^sea = from the beginning. 

Apart from the chiral zero modes, all eigenvalues of D come in complex conjugate 
pairs ibzAfe that can be ordered according to 

< Ai < As < . . . (2.4) 

The expectation values of these eigenvalues scale proportionally to 1/N. In partic- 
ular, the spectral densities at fixed topology, 

PkA^)= lim {6iz-XkN))^, 0<z<oo, (2.5) 

JV^oo 

as well as the corresponding joint distributions of the first n eigenvalues are well 

defined and analytically calculable [8,9]. Using these expressions, the expectation 
values of the first few eigenvalues and of their products can be obtained exactly, 
although in practice the formulae rapidly become so complicated that numerical 
methods are required. 
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3. Numerical simulation 



We assume that the reader is famihar with the standard formulations of lattice QCD 
and the Ginsparg- Wilson relation as discussed in refs. [16-20], for example. In this 
section we summarize the parameter choices that we have made and discuss what 
precisely is being compared with random matrix theory. 

3.1 Lattice theory 

We consider four-dimensional lattices of size L in all dimensions and impose periodic 
boundary conditions on the fields. The action of the SU(3) gauge field is taken to 
be the standard plaquette action with bare coupling go. As already mentioned, the 
lattice Dirac operator D should satisfy the Ginsparg-Wilson relation, 

j^D + L>75 = aDj^D, (3.1) 

in addition to the usual requirements such as locality. For this study we decided to 

use the Neuberger-Dirac operator [18], with shift parameter s = 0.4 [20], in view of 
its relative simplicity. The parameter a in eq. (3.1) is then given by 

a = a/{l + s), (3.2) 

where a denotes the lattice spacing (our notational conventions are as in refs. [10,20]). 

As usual we define the index i/ of D to be the difference n+ — n_ of the numbers 
of exact zero modes with positive and negative chirality. A well-known implication 
of the Ginsparg-Wilson relation is that u = Q, where 

Q = a^^g(x), g(x) = -iatr{75L)(x,x)}, (3.3) 

X 

can be taken as the definition of the topological charge of the gauge field [17]. In 
particular, the charge density q{x) is a local gauge-invariant expression in the link 
variables [20]. 

3.2 Numerical techniques 

At large volumes the low-lying non-zero eigenvalues of D are orders of magnitude 
smaller than the largest eigenvalues. It is hence important to choose a numerical 
implementation of the Neuberger-Dirac operator, where the approximation errors 
can be guaranteed to be sufficiently small for the spectrum to be obtained to the 
required level of precision. 
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Table 1. Simulation parameters 



lattice 




L/a 


ro/a 


L [fm] 


'meas 


Ai 


6.0 


12 


5.368 


1.12 


2452 


A2 


6.1791 


16 


7.158 


1.12 


1138 


Bo 


5.8458 


12 


4.026 


1.49 


2918 


Bi 


6.0 


16 


5.368 


1.49 


1001 


B2 


6.1366 


20 


6.710 


1.49 


963 


Co 


5.8784 


16 


4.294 


1.86 


1109 


Ci 


6.0 


20 


5.368 


1.86 


931 



The minmax polynomial approximation that was introduced in ref. [10] provides 
a solution to this problem. We also apply some of the other numerical techniques 
mentioned there. In particular, the eigenvalues of D are obtained by minimization 
of the Ritz functional [25,26] of the hermitian operators 

= P±DP±, P± = 1(1 ± 75), (3.4) 

in the positive and negative chirality sectors. 

3. 3 Simulation parameters 

We have simulated altogether 7 lattices with various sizes and lattice spacings (see 
table 1). There are three groups of lattices, labelled A, B and C, where the lattice 
size L is kept fixed in physical units. Within each of these groups, the bare coupling 
/? = 6/^0 is thus the only parameter that varies, which allows us to obtain a direct 
check on the dependence of the calculated observables on the lattice spacing. For 
the conversion to physical units we use the recent parametrization of the Sommer 
scale ro [27] by the ALPHA collaboration (cq. (2.6) of [28]) and set tq = 0.5 fm. 

In all these simulations the generation of the gauge-field configurations consumes 
a negligible amount of computer time. We have therefore performed many update 
cycles between subsequent configurations (typically 500 iterations of 1 heatbath and 
at least 6 over-relaxation updates of all link variables) so that they can be assumed 
to be statistically independent. The number of gauge fields in each ensemble quoted 
in the last column of table 1 coincides with the number of ""measurements'" of the 
topological charge and the lowest few eigenvalues of the Dirac operator. 
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3-4 Matching with random matrix theory 

The eigenvalues 7 of D lie on a circle in the complex plane, 

D^ = j^, 7=^(1-6^'^), (3.5) 

and come in complex conjugate pairs (if Im7 does not vanish). We cannot directly 
compare these eigenvalues with the eigenvalues of the random matrix T) since the 
latter lie on the imaginary axis. The important point to note is, however, that the 
radius of the circle diverges in the continuum limit and that the real parts of the 
eigenvalues 7 with I7I ^ 1/a rapidly go to zero in this limit. In practice we set (for 
these eigenvalues) 

A=|7| = ^{2(l-cos0)}^/' (3.6) 

and compare the distributions of the scaled eigenvalues 

z = ASy, V = L^ (3.7) 

with those of the scaled eigenvalues z = \N in the matrix model. Evidently the 
definition (3.6) is arbitrary to some extent, but other suggested formulae differ by 
terms of order aPX^ and are thus asymptotically equivalent. In particular, for the 
analysis of our data this ambiguity turns out to be numerically insignificant. 

Although wc use the same symbol for the proportionality constant E in eq. (3.7) as 
for the quark condensate in full QCD, it should be emphasized that in the present 
context S is just a free parameter with no obvious physical interpretation. The 
quark condensate is actually not even well defined in quenched QCD, since it is not 
possible to pass to the chiral limit in this theory. We thus prefer to regard S as an 
effective parameter that is to be determined, on any given lattice, by comparing the 
spectral distributions of D with those of D. 



4. Distribution of the topological charge 

We now first discuss the probability Pi, to find a gauge field with topological charge 
Q = I/. Random matrix theory does not provide any information on P^, and we 
are thus uniquely concerned with properties of QCD in this section. Basically we 
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wish to find out whether the charge distribution scales in the expected way as a 
function of the volume and the lattice spacing. If this is the case, it will then be 
more plausible that the eigenvalues of the Dirac operator behave coherently in the 
different sectors (as suggested by random matrix theory). 

4-1 Large volume limit 

The fact that the topological charge Q is given in terms of a local density q{x) allows 
us to derive an asymptotic formula for at large volumes. First note that since Q 
is integer- valued, we havef 

P^= r ^e-^^'-e-^W, F(0) = -ln{(e^^«)}. (4.1) 

Using the moment-cumulant transformation, the free energy F{9) may then be ex- 
panded in a series 

m = V^{-ir+'j^Cn, (4.2) 
Cn = a^'^-^ Yl • • • <l{^2n-i)q{0)r'' . (4.3) 

Xl,...,X2n-l 

In quenched QCD the connected correlation functions in this formula are evaluated 
in the pure gauge theory, where all particles are fairly heavy (the mass of the lightest 
glueball is around 1.6 GeV). As a result the cumulants C„ approach their infinite- 
volume limit exponentially fast, i.e. as soon as L is larger than 1 fm or so, they can 
be expected to be practically independent of the volume. 

If we now insert eq. (4.2) in eq. (4.1), it is clear that the integral is dominated by 
the saddle point at = in the large-volume limit. The asymptotic formula 



" {l + 0(y-i)} (4.4) 



is thus obtained, where = (Q^)- On the lattices Ai, . . . , Ci, the observed charge 

distributions are in fact statistically consistent with the leading term in eq. (4.4) (see 
fig. 1). From this point of view the lattices are hence in the large- volume regime. 



f As in the matrix model, expectation values in QCD at fixed charge Q = v are denoted by (. . .)^. 
Brackets without lower index imply an unconstrained expectation value and {(j)i{xi) . . . </'n(a;„))™" 
stands for the connected part of the full n-point correlation function of the fields <p\, . . . ,4>n- 
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Fig. 1. Histogram of the topological charge on lattice B2. The curve represents the 
leading term in the large-volume formula (4.4) with = (Q^) taken from table 2. 

4-2 Determination of the topological susceptibility 

The width of the charge distribution (or, cquivalently, the susceptibility x = (Q^) /^) 
may be calculated straightforwardly by averaging over the ensemble of gauge con- 
figurations that was generated. However, this procedure may not be safe in general, 
because the tails of the charge distribution are poorly sampled. 

To understand what the problem is, first note that the number Ni, of configura- 
tions with charge Q = v m a sample of A'^ configurations is distributed according to 
Poisson statistics. In particular, the variance of A^j^ is equal to its mean value NP^. 
If Ni, is small, the empirical probability N^i/N may consequently be quite different 
from the true probability P^- We thus compute (Q^) as a sum of two contributions, 



where fmax is the maximal charge such that N^, + N^^ is at least 10 for all v in the 
first sum. To evaluate the second sum, we insert the large-volume expression (4.4) 
for the exact distribution Pjy, setting to the naive estimate of (Q^), where any 
exceptional configurations with absolute charges \Q\ significantly larger than i^max 




(4.5) 
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Table 2. Topological susceptibility 





Ai 


A2 


Bo 


Bi 


B2 


Co 


Ci 


4 

Xro 


1.63(5) 
0.065(2) 


1.59(8) 
0.064(3) 


5.6(2) 
0.071(2) 


5.6(3) 
0.071(4) 


4.8(2) 
0.061(3) 


15.0(7) 
0.078(4) 


12.8(9) 
0.066(5) 



The errors quoted in the second row do not include the error on ro/a 



should be discarded (there were none in our samples). The second sum is actually 
a small correction to the first, and other estimates of could therefore be used at 
this point with little effect on the results listed in table 2.f 

The fact that x'^o comes out to be the same within errors on the lattices Ai,Bi 
and Ci (which have exactly the same lattice spacing but different volumes) is in line 
with the above argumentation that finite-volume effects in connected correlation 
functions of local fields should be small in the pure gauge theory when L > 1 fm. As 
far as the volume dependence is concerned, the charge distribution on the lattices 
that we have simulated thus behaves entirely according to expectations. 

4-3 Continuum limit 

It is often taken for granted that the division of the space of gauge fields into topo- 
logical sectors and the topological susceptibility have a well-defined meaning in the 
continuum limit of QCD. Beyond the semi-classical approximation, where the func- 
tional integral is expanded about the instanton solutions, the question remains un- 
decided, however, and is in fact difficult to pose in precise terms without reference 
to a regularization of the theory. In particular, the susceptibility cannot simply be 
defined as the integral over the two-point function {q{x)q{0)) of the charge density, 
because there is a non-integrable singularity at x = 0. 

In lattice QCD the space of fields is connected and the assignment of a topological 
charge to every lattice gauge field is consequently not unique. The definition (3.3), 
for example, depends on the choice of D. Close to the continuum limit, the integra- 
tion measure in the functional integral may, however, be increasingly supported on 
fields where the charge assignment is unambiguous. The topological sectors would 

f To the extent they can be compared, the results obtained in a similar study of the topological 
susceptibility by Del Debbio and Pica [29] , which appeared during the completion of the present 
paper, are compatible with those reported here. 
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then arise dynamically and any differences in the definition of the charge would be 
ultimately irrelevant. 

This picture suggests that the susceptibility should behave like a physical quantity 
of dimension 4, i.e. that the combination xfo should be independent of the lattice 
spacing, up to corrections of order a^. The data listed in the second row of table 2 
may actually be fitted by the linear expression x^o — cq + cio^, and our results 
are, therefore, statistically consistent with the existence of a well-defined continuum 
limit of the topological susceptibility. If we fit the data from the B lattices only, 
taking the error on ro/a into account [28], the value x^o ~ 0.059(5) is obtained at 
a = 0, but this number should evidently be used with caution, as it is determined 
by linear extrapolation in of only three data points. 



5. Comparison with random matrix theory 

As discussed in subsect. 3.4, the matching of the QCD spectra with random matrix 
theory involves an unknown scale S that must be determined from the data. This 
complication can be avoided by considering ratios of expectation values, and we now 
look at some of these first. 

5.1 Scale-independent tests 

In table 3 we list our results for the ratios {\k)i>/ {^j)^ on the lattices Ai, . . . , Ci. 
There is a significant cancellation of statistical errors when the ratios are formed, 
particularly for the larger values of k and j. For this reason we show the data for all 

1 < J < A; < 4, even though some of them arc related to each other. The figures in 
the last column of table 3 are the values of the ratios predicted by random matrix 
theory [8,9]. 

A number of observations can be made at this point: 

(1) For fixed physical volume (i.e. within the groups of lattices A, B or C), the data 
on each line are constant within 1 or, in rare cases, at most 2 standard deviations. 
In other words, there is no statistically significant dependence on the lattice spacing 
in these ratios. 

(2) In all topological sectors the results obtained on the lattices B agree with random 
matrix theory. There arc deviations at the level of 2 or 3 standard deviations in a 
few places, but this has to be so by the laws of statistics. 
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Table 3. Simulation results for the ratios {Xk)u / (Aj)i/ 







Ai 


A 

A2 


Bo 


Bi 


B2 


Co 


Ci 


KMi 





2/1 


2.29(4) 


2.28(6) 


2.71(6) 


2.73(10) 


2.56(10) 


2.77(12) 


3.01(14) 


2.70 





3/1 


3.29(7) 


3.25(9) 


4.45(11) 


4.59(18) 


4.43(18) 


4.65(22) 


4.89(23) 


4.46 





4/1 


4.07(8) 


4.00(11) 


6.25(15) 


6.55(26) 


6.15(26) 


6.67(31) 


6.99(34) 


6.22 





3/2 


1.44(1) 


1.42(2) 


1.65(2) 


1.68(3) 


1.73(3) 


1.68(4) 


1.63(4) 


1.65 





4/2 


1.78(2) 


1.75(2) 


2.31(3) 


2.40(5) 


2.40(5) 


2.41(6) 


2.32(6) 


2.30 





4/3 


1.24(1) 


1.23(1) 


1.40(1) 


1.43(2) 


1.39(2) 


1.44(2) 


1.43(2) 


1.40 


1 


2/1 


1.78(2) 


1.73(3) 


2.04(3) 


2.04(4) 


2.12(4) 


2.03(5) 


1.95(5) 


2.02 


1 


3/1 


2.35(3) 


2.27(4) 


3.08(4) 


3.08(6) 


3.23(7) 


3.13(8) 


2.97(7) 


3.03 


1 


4/1 


2.80(3) 


2.69(4) 


4.06(6) 


4.08(9) 


4.30(9) 


4.27(11) 


4.15(10) 


4.04 


1 


3/2 


1.32(1) 


1.31(1) 


1.51(1) 


1.51(2) 


1.53(2) 


1.55(2) 


1.53(2) 


1.50 


1 


4/2 


1.57(1) 


1.56(1) 


1.99(2) 


2.00(3) 


2.03(3) 


2.11(3) 


2.13(4) 


2.00 


1 


4/3 


1.19(1) 


1.19(1) 


1.32(1) 


1.32(1) 


1.33(1) 


1.36(2) 


1.40(2) 


1.33 


2 


2/1 


1.55(2) 


1.55(3) 


1.80(2) 


1.83(4) 


1.80(4) 


1.84(4) 


1.82(3) 


1.76 


2 


3/1 


1.98(3) 


1.95(4) 


2.55(3) 


2.60(6) 


2.51(6) 


2.66(6) 


2.66(5) 


2.50 


2 


4/1 


2.28(3) 


2.27(5) 


3.26(4) 


3.30(8) 


3.23(8) 


3.50(8) 


3.51(7) 


3.24 


2 


3/2 


1.28(1) 


1.26(1) 


1.42(1) 


1.42(2) 


1.39(2) 


1.44(2) 


1.46(2) 


1.42 


2 


4/2 


1.48(1) 


1.47(2) 


1.82(2) 


1.80(3) 


1.79(3) 


1.90(3) 


1.93(3) 


1.83 


2 


4/3 


1.16(1) 


1.16(1) 


1.28(1) 


1.27(1) 


1.29(1) 


1.32(1) 


1.32(1) 


1.29 



(3) The results on the lattices C are also matched by random matrix theory, although 
in this case some differences of up to 4 standard deviations are observed at = 1 and 
V = 2. It is still possible that these are caused by an unlikely statistical fluctuation, 
particularly so since they do not appear to follow any obvious systematic trend. 

(4) In the case of the lattices A, on the other hand, the data are in clear disagreement 
with random matrix theory. Since the ratios are practically independent of the lattice 
spacing, we conclude that we are seeing the onset of finite-volume effects that are 
specific to QCD and are not accounted for by random matrix theory. Such a regime 
must in any case eventually be reached when L is set to values below 1 fm. 

For illustration the ratios calculated on lattice B2 are plotted in fig. 2 together 
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Fig. 2. Comparison of simulation results for {\k)v/ {^j)v from lattice B2 (diamonds) 
with random matrix theory (horizontal bars) in the sectors with topological charge 
i/ = 0,l,2. 

with the values predicted by random matrix theory. This shows rather clearly that 
the observed matching is non-trivial and at a high level of precision. Evidently ex- 
pectation values of products of eigenvalues could also be considered, but it is our 
experience that ratios involving these are obtained with relatively large statistical er- 
rors, and comparisons with random matrix theory are consequently less compelling. 

5.2 Calculation ofT, 

In the present context S is considered to be a parameter with no independent phys- 
ical meaning whose value may be fixed by imposing a suitable normalization condi- 
tion. We may, for example, require the relation 

^X,)^CB^V = (z,)r^ (5.1) 

to hold exactly for some k and u, where Zk denotes the A;th scaled eigenvalue in the 
random matrix model (cf. subsect. 3.4). The numbers obtained in this way are listed 
in table 4. 

In the case of the lattices Ai and A2, we already know that the eigenvalue distri- 
butions in QCD are poorly matched by those in random matrix theory. The results 
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Table 4. Values of Ero determined from (Afc)^ 



V 


A; 


Ai 


A2 


Bo 


Bi 


B2 


Co 


Ci 


n 


1 

J. 






n 991 




U.ZfOUl ±o J 









2 


0.254(3) 


0.287(4) 


0.220(3) 


0.271(7) 


0.297(8) 


0.257(6) 


0.298(9) 





3 


0.291(2) 


0.332(3) 


0.221(2) 


0.266(5) 


0.283(5) 


0.253(5) 


0.302(7) 





4 


0.329(2) 


0.377(2) 


0.220(2) 


0.260(3) 


0.284(4) 


0.246(4) 


0.295(6) 


1 


1 


0.244(3) 


0.263(5) 


0.223(3) 


0.268(7) 


0.306(7) 


0.271(9) 


0.290(8) 


1 


2 


0.276(2) 


0.307(3) 


0.220(2) 


0.266(4) 


0.292(4) 


0.270(4) 


0.300(6) 


1 


3 


0.315(2) 


0.352(2) 


0.219(2) 


0.264(3) 


0.287(3) 


0.262(4) 


0.296(5) 


1 


4 


0.352(2) 


0.395(2) 


0.222(2) 


0.266(3) 


0.288(3) 


0.257(3) 


0.282(4) 


2 


1 


0.265(4) 


0.298(7) 


0.228(3) 


0.273(7) 


0.291(8) 


0.275(7) 


0.312(7) 


2 


2 


0.301(3) 


0.340(4) 


0.224(2) 


0.262(5) 


0.284(5) 


0.262(6) 


0.301(5) 


2 


3 


0.335(2) 


0.382(3) 


0.223(2) 


0.262(4) 


0.289(4) 


0.258(4) 


0.293(5) 


2 


4 


0.374(2) 


0.425(3) 


0.226(2) 


0.267(3) 


0.291(3) 


0.253(3) 


0.287(4) 


The 


errors 


> quoted in 


this table do not include the error on 


ro/a 







for S consequently show a strong dependence on k and v, and any normalization 
convention that one may adopt is therefore rather arbitrary. The situation looks 
much better on the lattices Bo, . . . , Ci, where the calculated values of S are prac- 
tically the same for all k and u (i.e. in each column of table 4). Up to the quoted 
errors, S is thus consistently determined in these cases. 

There is, however, no reason to expect the dimensionless combination Erg to be 
independent of the lattice spacing and the lattice size. In particular, at fixed lattice 
spacing, Srg appears to be monotonically increasing with the volume. Close to the 
continuum limit, S presumably renormalizes like the scalar quark density S, so that 

ZsJ^r^ = hin/A, AL) + 0(0^), (5.2) 

where Zs{go,aiJ,) denotes a renormalization constant, the renormalization mass 
and A the associated renormalization-group invariant scale. The analogy with full 
QCD suggests this, and, more importantly, the fact that the higher-order cumulants 
of S (which renormalize multiplicatively) are directly related to the spectral density 
of the Dirac operator [2]. 
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Fig. 3. Probability distribution pi{t) of the local magnitude t = \'>pi(x)\'^ V of the 
first normalized eigenvector of P+D^ DP+ in the charge Q = sector. The solid and 
dotted lines are the distributions obtained on the lattices Ci and B2 respectively, while 
the dashed curve is what would be expected for a random vector with unit norm. 

The renormalization factor Zs connecting the lattice theory to the MS scheme 
of dimensional regularization has been calculated to one-loop order of perturbation 
theory [30,31]. This formula is unfortunately not very useful in the present context, 
because it largely underestimates the true value of the renormalization constant 



which is obtained from non-perturbativc scaling studies [32]. As a consequence, and 
since Zs is (for s = 0.4) currently only known at this value of /3, we are unable to 
study the scaling behaviour of Z^Srg. More extensive simulations would in any case 
be required for a reliable extrapolation to the continuum limit (if eq. (5.2) holds). 

5. 3 Local magnitude of the low modes 

In correlation functions of local operators such as the pseudo-scalar and the scalar 
quark densities, the contribution of the low modes of the Dirac operator can be 
large, owing to the presence of uncancelled factors of 1/ Afc . The actual size of these 



Zs 



yu=2GeV ~ 



1.43(11) at P = 6.0, 



(5.3) 
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contributions also depends on the magnitude \ijjk{x)\^ of the associated normahzed 
eigenfunctions ipk{x) at the points x where the operators sit. 

In fig. 3 we show tlie distribution of the local magnitude of the lowest mode in the 
vacuum sector. Although an average over many gauge configurations is plotted, this 
curve is actually fairly universal, i.e. fluctuations are small and very similar results 
are obtained for the higher modes and in the other topological sectors. Random ma- 
trix theory does not refer to an underlying local structure and is hence not expected 
to be relevant here. In the two chirality sectors, the eigenvectors of WD are in fact 
uniformly distributed random vectors, and the associated distribution (which is also 
displayed in the figure) is clearly different from the one in QCD. 

An interesting aspect of the QCD distribution shown in fig. 3 is that it decreases 
only relatively slowly at large local magnitudes t. The inset in the figure shows 
this in greater detail. In particular, the probability to find a point on the lattice 
where t > 3 is more than 5%, and there is still a non- negligible probability of nearly 
0.6% for having t > 10. Correlation functions of local operators in the e-regime may 
consequently suffer from large statistical fluctuations that derive from the presence 
of exceptionally low eigenvalues of the Dirac operator in combination with accidental 
"bumps" in the associated wave functions at the positions of the operators. 



6. Conclusions 



The numerical simulations reported in this paper lend further support to the propo- 
sition that the low-lying eigenvalues of the Dirac operator in the e-regime of QCD are 
distributed according to the chiral unitary random matrix model. Detailed agree- 
ment has been observed, in different topological sectors, on all lattices with linear 
extent L larger than about 1.5 fm. 

Random matrix behaviour sets in rather rapidly when going from small to large 
volumes. To some extent this can be understood by noting that the average of the 
lowest eigenvalue in the vacuum sector is inversely proportional to the volume. On 
the lattices Ai, Bi and Ci, for example, the renormalized spectral gap 

A = Zgi(Ai)o (6.1) 

in the MS scheme at /i = 2GeV is approximately equal to 91, 23 and 8MeV respec- 
tively. Clearly the e-regime has been safely reached on Ci, while this is not so in 
the case of the lattice Ai. 
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Our results also suggest that the lattice theory scales coherently to the continuum 
limit in all topological sectors (as is usually assumed). More extensive simulations 
would, however, be needed, to be able to take the continuum limit of the topological 
susceptibility and other quantities with confidence. This is technically possible but 
requires significant computer resources, since the numerical effort in these calcula- 
tions tends to grow roughly proportionally to (L/a)^(L/ro)^. 

We are indebted to Poul Damgaard, Pilar Hernandez, Karl Jansen, Laurent Lel- 
louch, Ferenc Niedermayer and Giancarlo Rossi for helpful discussions. The simu- 
lations were performed on PC clusters at DESY-Hamburg, the Institut fiir Theo- 
retische Physik at the University of Bern, the Max-Planck-Institut fiir Physik in 
Munich, the Max-Planck-Institut fiir Plasmaphysik in Garching, and the Leibniz- 
Rechenzentrum der Bayerischen Akademie der Wissenschaften. We wish to thank all 
these institutions for supporting our project and the staff of their computer centres 
(particularly Peter Breitenlohner, Isabel Campos and Andreas Gellrich) for technical 
help. L. G. was supported in part by the EU under contract HPRN-CT-2000-00145 
Hadrons/LatticeQCD. 
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